1. Loading data

1.1. phyloseq obejct

1.2. qPCR data (controls)

2. QC

QC1. How many samples failed sequencing

QC2. How were changes in read stats and host DNA proportion?

QC3. How were the extraction controls

QC4. Prevalence / abundance filtering - red flag

3. Analysis

A0. Calculation of alpha-diversity indices

A1. Host DNA, bacterial DNA and % host

A2. Modeling of sequencing results

A3. Taxa alpha diversity

A4. Taxa beta diversity

Intermediate results

A5. DA analysis for taxa

A6. Decontam

A7. LM of function alpha diversity (BPI)

A8. permanova of function alpha diversity

A9. DA for function

Data inputs

Meta data

• qPCR - bacteria

• qPCR - human

• qPCR host %

• Raw reads

• final reads

• sequencing host %

• library prep failure status

• Raw reads

• subject_id

• treatment

• sample_type

• subject_id

Sequencing result

• samples

• controls

QC table by treated (binary)

Changes in matrices were observed

#sequencing result by sample type and control (1/0)
options(dplyr.summarise.inform = FALSE)

sample_data %>% data.frame() %>% 
        group_by(sample_type, treated) %>% 
        summarise(N = n(),
                  `Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2),nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
        ) %>%
        rename(`Sample type` = sample_type, Treated = treated) %>%
        data.frame(check.names = F) %>% mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")

QC table by treatment methods

Changes were sample type * treatment specific

sample_data %>% data.frame() %>% 
        #dplyr::filter(sample_type %in% c("Sputum", "nasal_swab", "BAL")) %>% 
        group_by (sample_type, treatment) %>%
        summarise(N = n(),
              `Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
              `Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
              `Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
              `Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
        ) %>% data.frame(check.names = F) %>% 
        arrange(sample_type, treatment) %>%
        rename(`Sample type` = sample_type, Treatment = treatment) %>%
        mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")

Figure of reads by treatment (z-score)

Changes were sample type * treatment specific

# Summary figures - facet and z-score -------------------------------------

sample_data %>% 
        subset(., !is.na(.$subject_id)) %>%
        data.frame() %>%
        gather(feature, value, Raw_reads:sequencing_host_prop) %>%
        group_by(feature, sample_type) %>% 
        subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop")) %>%
        mutate(z_score = scale(value),
               feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
        ggplot(aes(x = treatment, y = z_score, fill = treatment)) +
                geom_boxplot(lwd = 0.2) +
                guides(fill=guide_legend(title="Treatment", nrow = 1)) +
                facet_grid(sample_type~feature) +
                xlab("Treatment") +
                ylab("Z score") +
                theme_classic(base_family = "serif", base_size = 14) +
                guides( x =  guide_axis(angle = 90)) + 
                theme(legend.position = "top") +
                scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment") #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6

Results:

2.1. There were no differences in raw reads.

2.2. However, final reads increased after some treatment, and host DNA proportion decreased

Species richness of controls

Some possible contaminants were identified in extraction controls

#Loading theoretical mock community
zymo_mock <- read_excel("/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_sicas2/data_raw/DAR_20210929_zymo_mock_data.xlsx") %>%
        data.frame(row.names = T) %>% rename(mock_theoretical = Mock) %>% mutate(mock_theoretical = mock_theoretical/100) %>%
        merge_phyloseq(otu_table(., taxa_are_rows = T), tax_table(phyloseq$phyloseq_count))

phyloseq_mock <- rbind(c("mock_theoretical", "Mock theoretical", "-")) %>% data.frame() %>%
        column_to_rownames(var = "X1") %>% rename(sample_type = X2, treatment = X3) %>% #making sample_data of mock community
        merge_phyloseq(sample_data(.), zymo_mock)

phyloseq_control_rel <- subset_samples(phyloseq$phyloseq_rel, sample_type == "Mock" | sample_type == "Neg.") #adding data of controls
sample_data(phyloseq_control_rel)$treatment <- sample_data(phyloseq_control_rel)$treatment %>% as.character()
sample_data(phyloseq_control_rel)$sample_type <- sample_data(phyloseq_control_rel)$sample_type %>% as.character()
phyloseq_control_rel <- merge_phyloseq(phyloseq_control_rel, phyloseq_mock) 


#Species richness of each control groups
sample_data(phyloseq_control_rel)$S.obs <- rowSums(t(otu_table(phyloseq_control_rel)) != 0)
sample_data(phyloseq_control_rel)$sample_type <- 
        factor(sample_data(phyloseq_control_rel)$sample_type, levels = c("Mock theoretical", "Mock", "Neg."))
sample_data(phyloseq_control_rel)$teratment <-
        factor(sample_data(phyloseq_control_rel)$treatment, levels = c("-", "Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"))
phyloseq_control_rel %>%
        sample_data() %>%
        #mutate(sample_type = factor(sample_type, levels = c("Mock", "Neg.")),
        #       treatment = factor(treatment, levels = c("Theoretical", "Untreated", "Benzonase", "Host zero", "Molysis", "QIAamp"))) %>%
        group_by(sample_type, treatment) %>%
        summarise(Mean = mean(S.obs),
                  SD = sd(S.obs),) %>%
        kbl(format = "html", caption = "Species richness of controls") %>% 
        kable_styling(full_width = 0, html_font = "serif")

Bar plot of controls

Some possible contaminants were identified in extraction controls

Some changes visible at postive control….

#Manipulating phyloseq - only top 10 
tax_table(phyloseq_control_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Bar plot of control data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ sample_type, scales= "free_x", nrow=1)

#there could be opportunistic pathogens...

Bar plot of controls (Positive)

Some possible contaminants were identified in extraction controls

Gram negatives were fragile to depletion method at postivive control

#Manipulating phyloseq - only top 10 

phyloseq_control_rel %>% 
        subset_samples(., sample_type == "Mock") %>% 
        tax_table() %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Postive controls") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

#there could be opportunistic pathogens...

Issue 1 - gram negative

Benzonse and Host zero depleted all the gram negative strains

Others decrased gram negatives a lot, but not became zero

# gram stain data


phyloseq_control_rel %>% 
        subset_samples(., sample_type == "Mock") %>% 
        tax_table() %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="Gram") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain in Zymo mock") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Gram-stain")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

#there could be opportunistic pathogens...
                                
#Manipulating phyloseq - only top 10 

sample_data(phyloseq_control_rel) <- cbind(phyloseq_control_rel %>%
              sample_data %>%
              data.frame(),
      phyloseq_control_rel %>% 
        otu_table %>%
        data.frame %>%
        subset(., rownames(.) %in% head(taxa_sums(subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)) %>%
        t()
)


sample_data(phyloseq_control_rel) %>%
        data.frame() %>%
        subset(., !is.na(.$Escherichia_coli)) %>% 
        group_by(sample_type, treatment) %>% 
        summarise(N = n(),
                  `<i>Escherichia coli</i><br>(median [IQR])` = paste(format(round(median(Escherichia_coli),4),nsmall = 4, big.mark = ","), " [", format(round(quantile(Escherichia_coli, 0.25),4), nsmall = 4, big.mark = ","), ", ", format(round(quantile(Escherichia_coli, 0.75),4), nsmall = 4, big.mark = ","), "]", sep = ""),
                  `<i>Pseudomonas aeruginosa</i><br>(median [IQR])` = paste(format(round(median(Pseudomonas_aeruginosa_group),4),nsmall = 4, big.mark = ","), " [", format(round(quantile(Pseudomonas_aeruginosa_group, 0.25),4), nsmall = 4, big.mark = ","), ", ", format(round(quantile(Pseudomonas_aeruginosa_group, 0.75),4), nsmall = 4, big.mark = ","), "]", sep = ""),
                  `<i>Salmonella enterica</i><br>(median [IQR])` = paste(format(round(median(Salmonella_enterica),4),nsmall = 4, big.mark = ","), " [", format(round(quantile(Salmonella_enterica, 0.25),4), nsmall = 4, big.mark = ","), ", ", format(round(quantile(Salmonella_enterica, 0.75),4), nsmall = 4, big.mark = ","), "]", sep = "")
        ) %>%
        rename(`Sample type` = sample_type) %>%
        data.frame(check.names = F) %>% mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")

Issue 2 - Positive controls contaminants

Some possible contaminants were identified in most of samples

This could be 1) background contamination or 2) cross-contamination from kingfisher. Most of these are gram positives. Negative controls should be double-checked

#Manipulating phyloseq - only top 10 
phyloseq_control_rel_contam <- subset_taxa(phyloseq_control_rel , !(taxa_names(phyloseq_control_rel) %in% head(taxa_sums(subset_samples(phyloseq_control_rel,
                                                                                         sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10))
)

phyloseq_control_rel_contam <- subset_taxa(phyloseq_control_rel_contam, taxa_sums(phyloseq_control_rel_contam) != 0)
phyloseq_control_rel_contam <- subset_samples(phyloseq_control_rel_contam, sample_type != "Neg." & S.obs != 0)


tax_table(phyloseq_control_rel_contam) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel_contam) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel_contam
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Contaminants in Zymo mock") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

Negative controls

Contaminants of positive and negative control do not match

It seems without host DNA, gram-negatives are vulnerable to depletion methods.

These negative contaminants highly likely introduced after-depletion

#Manipulating phyloseq - only top 10 

tax_table(phyloseq_control_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "Neg.") %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Barplot of neg. data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

phyloseq_control_rel %>%
        subset_samples(., sample_type == "Neg.") %>%
        plot_bar(., fill="Gram") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram-stain of negative data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

Samples - gram-stain?

Nasal swab had low gram positive/pattern consisted after depletion

BAL showed Similar gram - / + ratio

Sputum showed high decrease in gram negative bacteria

Freeze/thaw cycle could be associated

Currently no further analysis is possible

#Manipulating phyloseq - only top 10 

tax_table(phyloseq$phyloseq_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq$phyloseq_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq$phyloseq_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "Nasal") %>%
        plot_bar(., fill="Gram") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain of nasal samples") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

tax_table(phyloseq$phyloseq_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq$phyloseq_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq$phyloseq_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "BAL") %>%
        plot_bar(., fill="Gram") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain of BAL samples") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

tax_table(phyloseq$phyloseq_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq$phyloseq_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq$phyloseq_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "Sputum") %>%
        plot_bar(., fill="Gram") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain of sputum samples") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

Histogram of prelanence taxa

No prevalence or abundance filtering (each experimental group is 5% of total sample)

#Calculation of sample prevalence, standard deviation, median abundance across all samples for all bugs and making into a table.
#
#•  In initial analysis we will not perform prevalence or abundance filtering (though we may consider this for secondary differential abundance analyses to manage p (features) > n (sample size) problem and issues with multiple hypothesis correction)
taxa_qc <- data.frame("species" =  otu_table(phyloseq$phyloseq_rel) %>% t() %>% colnames(),
        "prevalence" = ifelse(phyloseq$phyloseq_count %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
        "mean_rel_abd" = phyloseq$phyloseq_rel %>% otu_table() %>% t() %>% colMeans(na.rm = T) #mean relativ abundacne 
)

hist(log10(taxa_qc$prevalence), xlab = "log10(Taxa prevalence)", main = "Histogram of prevalence of taxa")

Histogram of mean abundance

hist(log10(taxa_qc$mean_rel_abd), xlab = "log10(Mean relative abundance)", main = "Histogram of mean relative abundance")

Red flag taxa

Taxa with low prevalences were red-flagged

red_flag_taxa <- data.frame(species = taxa_qc$species,
                          red_flag_prev_abd = ifelse(taxa_qc$prevalence < otu_table(phyloseq$phyloseq_rel) %>% t %>% rownames() %>% length * 0.05 & taxa_qc$mean_rel_abd < quantile(taxa_qc$mean_rel_abd, 0.75), 1, 0))
red_flag_taxa

qPCR result

#2A: Change in total DNA (qPCR)
f2a <- ggplot(sample_data, aes(x = sample_type, y = log10(DNA_host_nondil + DNA_bac_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        ylab("log<sub>10</sub>(qPCR total DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif") + 
        labs(tag = "A") +
        #scale_x_discrete(label = c( "Mock", "Neg.", "BAL", "Nasal", "Sputum")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))


#2B: Change in human DNA (qPCR)
f2b <- ggplot(sample_data, aes(x = sample_type, y = log10(DNA_host_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) +
        ylab("log<sub>10</sub>(qPCR host DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif")+ 
        labs(tag = "B") +
        #scale_x_discrete(label = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#2C: Change in 16S DNA (qPCR)
f2c <- ggplot(sample_data, aes(x = sample_type, y = log10(DNA_bac_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) +
        ylab("log<sub>10</sub>(qPCR bacterial DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif")+ 
        labs(tag = "C") +
        #scale_x_discrete(label = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))

#2D. Change in % host (qPCR)
f2d <- ggplot(sample_data, aes(x = sample_type, y = host_proportion)) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) +
        ylab("Host DNA ratio") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif") + 
        labs(tag = "D") +
        #scale_x_discrete(label = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))

#output for markdown
ggarrange(f2a, f2b, f2c, f2d, common.legend = T , align = "hv")

Figure 2. qPCR result of host depletion study. A. Total DNA B. Host DNA C. Bacterial DNA D. Host %

Sequencing result

f3a <- ggplot(sample_data, aes(x = sample_type, y = log10(Raw_reads))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        theme_classic (base_size = 12, base_family = "serif") + 
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type") +
        theme(axis.title.y = element_markdown(),
              plot.tag = element_text(size = 15)) +
        ylab("log<sub>10</sub>(raw reads)") +
        labs(tag = "A") +
        guides(fill = guide_legend(nrow = 1))

#   - Host_mapped


f3b <- ggplot(sample_data, aes(x = sample_type, y = log10(Host_mapped))) +
        theme_classic (base_size = 12, base_family = "serif")+ 
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type")+
        theme(axis.title.y = element_markdown(),
              plot.tag = element_text(size = 15)) +
        ylab("log<sub>10</sub>(host reads)") +
        labs(tag = "B") +
        guides(fill = guide_legend(nrow = 1))


#   - % Host (we have used Host_mapped/Raw_reads in prior papers)

#   - Final_reads

f3c <- ggplot(sample_data, aes(x = sample_type, y = log10(Final_reads))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type")+
        ylab("log<sub>10</sub>(final reads)") +
        theme_classic (base_size = 12, base_family = "serif") + 
        theme(axis.title.y = element_markdown(),
              plot.tag = element_markdown(size = 15)) +
        labs(tag = "C") +
        guides(fill = guide_legend(nrow = 1))


#   - % Host (we have used Host_mapped/Raw_reads in prior papers)
f3d <- ggplot(sample_data, aes(x = sample_type, y = sequencing_host_prop)) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        theme_classic (base_size = 12, base_family = "serif")+ 
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type")+
        theme(axis.title.y = element_markdown(),
              plot.tag = element_text(size = 15)) +
        ylab("Host ratio by sequencing") +
        labs(tag = "D") +
        guides(fill = guide_legend(nrow = 1))



ggarrange(f3a, f3b, f3c, f3d, common.legend = T, align = "hv")

Figure 3. Sequencing result of host depletion study. A. Total DNA B. Host DNA C. Bacterial DNA D. Host %

Library failure - ANOVA

Some samples failed in library prep. What type of sample were fragile to treatments?

glm ( library fail ~ sample_type + treatment + sample_type * treatment + subject_id )

glmer(lib_failed ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        Anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>Chisq)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Library failure

glm ( sequencing fail ~ sample_type + treatment + sample_type * treatment + subject_id )

Nasals were fragile to lyPMA and Molysis

glmer(lib_failed ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Sequencing failure

Modeling of sequencing failure were not available due to low number of cases.

BAL079 - control & lyPMA failed sequencing.

sample_data(phyloseq$phyloseq_count) %>% data.frame %>% mutate(sequencing_fail = (S.obs == 0)) %>%
glmer(sequencing_fail ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = .) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

log10(Final reads) - ANOVA

Which methods was effective in increasing the final reads?

Interaction term was significant

lmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

log10(Final reads)

Which methods was effective in increasing the final reads?

lmer( log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id )

Except lyPMA, every methods increased final reads

lmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Host ratio ANOVA

Which methods was effective in lowering host %

Interaction term was significant

lmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Host ratio

Which methods was effective in lowering host %

lmer( Host DNA ratio vs sample_type + treatment + sample_type * treatment + (1|subject_id) )

Host zero was effect to to all types. Molysis was effective to Nasal and sputum. QIAamp was effective for Nasal only.

lmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>%
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x))  %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Gram negatives - ANOVA

Which methods was changed gram-strain ratio?

Square root transformation was required

Interaction term was significant

hist(sample_data %>% data.frame %>% .$gram_neg)

hist(sample_data %>% data.frame %>% .$gram_neg %>% sqrt())

lmer(sqrt(gram_neg) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Gram negatives

Which method biased gram positive-negative ratio

lmer( Gram-negative ratio vs sample_type + treatment + sample_type * treatment + (1|subject_id) )

Some treatment (commercial) changed gram negative proportion

lmer(sqrt(gram_neg) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>%
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x))  %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Species richness (all samples & interaction term) - ANOVA

Interaction term was significant

sample_data <- sample_data(phyloseq$phyloseq_count) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))

lmer_sob <- lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data)
lmer_sob %>% 
        anova() %>% 
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness (all samples & interaction term)

Increase at sputum was at every treatment Postive and negative control showed no changes

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness - stratified (Pos + Neg)

No treatment increased species richenss - after adjusting sequencing depth. With mock community except lyPMA, treatments showed they even reduced the possible contaminants. Need to observe alpha diversity of positive controls

lm(S.obs ~ sample_type * treatment + log10 (Final_reads), data = subset(sample_data, sample_data$sample_type == "Neg." | sample_data$sample_type == "Mock" )) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness - stratified (NS + Pos + Neg)

Molysis and host zero may incrased speciess richness of Nasal Data include nasal swab, positive depletion, and negative depletion

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal" | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness (BAL + Pos + Neg)

No changes observed

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "BAL"  | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness (sputum + Pos + Neg)

Benzonase may incrased speciess richness of sputum

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum"  | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inv Simp - ANOVA

Final reads did not affect inverse Simpson

lmer_invsimpson <- lmer(data_invsimpson ~ sample_type * treatment + log10(Final_reads) + (1|subject_id), data = sample_data)

lmer_invsimpson %>% 
        anova() %>% 
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

lmer_invsimpson <- lmer(data_invsimpson ~ sample_type * treatment + (1|subject_id), data = sample_data)

Simpson (all samples & interaction term)

Sputum after treatment showed differences - stratified analysis is required

#Simpson

lmer_invsimpson %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (Untreateds)

Inverse Simpson ~ sample_type + log10 (Final_reads) + (1|original_sample)

Mock community treated with lyPMA showed cahnges in alpha diveresity

lmer(data_invsimpson ~ sample_type * treatment + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (NS + Pos + Neg)

Inverse Simpson ~ sample_type + log10 (Final_reads) + (1|original_sample)

Mock community treated with lyPMA only showed changes in alpha diveresity.

lmer(data_invsimpson ~ sample_type * treatment + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal" | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (BAL + Pos + Neg)

Nothing changed in BAL

lmer(data_invsimpson ~ sample_type * treatment + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL" |sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (spt + Pos + Neg)

Sputum changed after some treatment - but their changes were not treatment global.

lmer(data_invsimpson ~ sample_type * treatment + log10(Final_reads) + + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum" |sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

PCoA based on Bray-Curtis (all samples)

Based on distances, it seems like some of the negative controls were affected by some samples. Meanwhile, Positive controls (mock community) were close to BAL samples.

ordinate(phyloseq_rel_nz,  method = "PCoA", distance = "bray") %>%
        plot_ordination(phyloseq_rel_nz, ., col = "treatment", shape = "sample_type" ) + 
        #scale_color_viridis(discrete = 6, name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) +
        scale_color_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_shape(name = "Sample type", labels = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        geom_point(size = 3) +
        theme_classic (base_size = 12, base_family = "serif") +
        theme(plot.tag = element_text(size = 15), legend.spacing = unit(0, 'cm'), legend.key.height = unit(0.4, "cm")) + #legend.position = c(0.9, 0.4)
        labs(tag = "E")

PCoA based on Jaccard (all samples)

Jaccard dissimilarities (presenece and absence) showed BAL and Mock communities are distant. Some samples may have some overlaps

ordinate(phyloseq_rel_nz,  method = "PCoA", distance = "jaccard") %>%
        plot_ordination(phyloseq_rel_nz, ., col = "treatment", shape = "sample_type" ) + 
        #scale_color_viridis(discrete = 6, name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) +
        scale_color_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_shape(name = "Sample type", labels = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        geom_point(size = 3) +
        theme_classic (base_size = 12, base_family = "serif") +
        theme(plot.tag = element_text(size = 15), legend.spacing = unit(0, 'cm'), legend.key.height = unit(0.4, "cm")) + #legend.position = c(0.9, 0.4)
        labs(tag = "E")

Stratified Bray beta diversity (Mock)

Some treatment made samples distant to the theoretical composition

ordinate(subset_samples(phyloseq_control_rel, sample_type != "Neg."), method = "PCoA", distance = "bray") %>%
        plot_ordination(phyloseq_control_rel, ., col = "treatment") +
        #scale_color_viridis(discrete = 6, name = "Treatment", labels = c("Mock theoretical", "Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) +
        scale_color_manual(values = c("black", "#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"),
                           name = "Treatment",
                           breaks = c("-", "Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"),
                           labels = c("Mock theoretical", "Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        #scale_shape(name = "Sample type", labels = c("Mock theoretical", "Mock")) +
        geom_point(size = 3) +
        theme_classic (base_size = 12, base_family = "serif") +
        theme(plot.tag = element_text(size = 15), legend.spacing = unit(0, 'cm'), legend.key.height = unit(0.4, "cm"))  #legend.position = c(0.9, 0.4)

        #labs(tag = "E")

Stratified Jaccard beta diversity (Mock)

Some treatment made samples distant to the theoretical composition

ordinate(subset_samples(phyloseq_control_rel, sample_type != "Neg."), method = "PCoA", distance = "jaccard") %>%
        plot_ordination(phyloseq_control_rel, ., col = "treatment") +
        scale_color_manual(values = c("black", "#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"),
                           name = "Treatment",
                           breaks = c("-", "Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"),
                           labels = c("Mock theoretical", "Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        #scale_shape(name = "Sample type", labels = c("Mock theoretical", "Mock")) +
        geom_point(size = 3) +
        theme_classic (base_size = 12, base_family = "serif") +
        theme(plot.tag = element_text(size = 15), legend.spacing = unit(0, 'cm'), legend.key.height = unit(0.4, "cm"))  #legend.position = c(0.9, 0.4)

        #labs(tag = "E")

Bar plot of mock community

Some bugs greatly decreased after some treatment

Differential abundance analysis should be conducted.

#Manipulating phyloseq - only top 10 
phyloseq_control_rel%>%
        subset_samples(sample_type != "Neg.") %>%
        tax_table(.) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel%>%
        subset_samples(sample_type != "Neg.") ) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 20)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel%>%
        subset_samples(sample_type != "Neg.") 
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Bar plot of positive controls") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 20 species")) +
        facet_wrap (~ treatment, scales= "free_x", nrow=1)

#there could be opportunistic pathogens...

Stratified Bray beta diversity (negative)

However, it seems there is no differences between negative controls..

ordinate(subset_samples(phyloseq_control_rel, sample_type == "Neg."), method = "PCoA", distance = "bray") %>%
        plot_ordination(phyloseq_control_rel, ., col = "treatment") +
        scale_color_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"),
                           name = "Treatment",
                           breaks = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"),
                           labels = c("Untreated","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        #scale_shape(name = "Sample type", labels = c("Mock theoretical", "Mock")) +
        geom_point(size = 3) +
        theme_classic (base_size = 12, base_family = "serif") +
        theme(plot.tag = element_text(size = 15), legend.spacing = unit(0, 'cm'), legend.key.height = unit(0.4, "cm"))  #legend.position = c(0.9, 0.4)

        #labs(tag = "E")

Beta diversity boxplot

Distances between samples within each subject. Mean distance between control <-> treatment for each subject

#distances of betadiversity - boxplots
bray_dist_long <- distance(phyloseq_rel_nz, method="bray") %>% as.matrix() %>% melt_dist() #making long data of distance matrices

#Adding sample type and treatment name. 
#this can be also done by merging metadata into the `bray_dist_long`
names <- data.frame(str_split_fixed(bray_dist_long$iso1, "_", 3))
names2 <- data.frame(str_split_fixed(bray_dist_long$iso2, "_", 3))
bray_dist_long$sample_id_1 <- paste(names$X1, names$X2, sep = "_")
bray_dist_long$method_1 <- ifelse(grepl("lyPMA", bray_dist_long$iso1),"lypma", 
                                         ifelse(grepl("ben", bray_dist_long$iso1),"benzonase", 
                                                ifelse(grepl("host", bray_dist_long$iso1),"host_zero", 
                                                       ifelse(grepl("qia", bray_dist_long$iso1),"qiaamp", 
                                                              ifelse(grepl("moly", bray_dist_long$iso1),"molysis", 
                                                                     "control")))))


#Adding data for iso 2 also should be done
bray_dist_long$sample_id_2 <- paste(names2$X1, names2$X2, sep = "_")
bray_dist_long$method_2 <-ifelse(grepl("lyPMA", bray_dist_long$iso2),"lypma", 
                                        ifelse(grepl("ben", bray_dist_long$iso2),"benzonase", 
                                               ifelse(grepl("host", bray_dist_long$iso2),"host_zero", 
                                                      ifelse(grepl("qia", bray_dist_long$iso2),"qiaamp", 
                                                             ifelse(grepl("moly", bray_dist_long$iso2),"molysis", 
                                                                    "control")))))


#subsetting distances of my interest
bray_dist_long$sample_id_1 <- ifelse(grepl("pos", bray_dist_long$sample_id_1, ignore.case = T),"Mock", 
                                 ifelse(grepl("neg|n_", bray_dist_long$sample_id_1, ignore.case = T),"Neg.",
                                        bray_dist_long$sample_id_1))
bray_dist_long$sample_id_2 <- ifelse(grepl("pos", bray_dist_long$sample_id_2, ignore.case = T),"Mock", 
                                 ifelse(grepl("neg|n_", bray_dist_long$sample_id_2, ignore.case = T),"Neg.",
                                        bray_dist_long$sample_id_2))


path_bray_dist_long_within_sampleid_from_control <- subset(bray_dist_long, bray_dist_long$sample_id_1 == bray_dist_long$sample_id_2) # data within samples

path_bray_dist_long_within_sampleid_from_control <- subset(path_bray_dist_long_within_sampleid_from_control,
                                                           path_bray_dist_long_within_sampleid_from_control$method_1 != path_bray_dist_long_within_sampleid_from_control$method_2) # remove irrelevant association

path_bray_dist_long_within_sampleid_from_control <- subset(path_bray_dist_long_within_sampleid_from_control,
                                                           path_bray_dist_long_within_sampleid_from_control$method_1 != path_bray_dist_long_within_sampleid_from_control$method_2) # remove irrelevant association

path_bray_dist_long_within_sampleid_from_control <- subset(path_bray_dist_long_within_sampleid_from_control, (path_bray_dist_long_within_sampleid_from_control$method_1 == "control") + (path_bray_dist_long_within_sampleid_from_control$method_2 == "control") != 0)

path_bray_dist_long_within_sampleid_from_control$treatment <- path_bray_dist_long_within_sampleid_from_control$method_1

path_bray_dist_long_within_sampleid_from_control$treatment <- ifelse(path_bray_dist_long_within_sampleid_from_control$treatment == "control", path_bray_dist_long_within_sampleid_from_control$method_2, path_bray_dist_long_within_sampleid_from_control$treatment) #Setting key method



path_bray_dist_long_within_sampleid_from_control$sample_type <- ifelse(grepl("NS", path_bray_dist_long_within_sampleid_from_control$iso1), "Nasal",
                                                                  ifelse(grepl("CFB", path_bray_dist_long_within_sampleid_from_control$iso1), "Sputum",
                                                                         ifelse(grepl("BAL", path_bray_dist_long_within_sampleid_from_control$iso1), "BAL",
                                                                                ifelse(grepl("pos|POS", path_bray_dist_long_within_sampleid_from_control$iso1), "Mock",
                                                                                       ifelse(grepl("neg|N_EXT", path_bray_dist_long_within_sampleid_from_control$iso1), "Neg.",NA)))))

path_bray_dist_long_within_sampleid_from_control %>%
        mutate(across(sample_type, factor, levels=c( "Neg.", "Mock", "BAL", "Nasal","Sputum"))) %>%
ggplot(aes(y = dist, fill = treatment)) +
        geom_boxplot() +
        #scale_fill_manual(values = c(viridis(6)[2:6])) +
        scale_fill_manual(values = c("#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        ylab("Sample pair distances") +
        theme_classic (base_size = 12, base_family = "serif") + 
        theme(plot.tag = element_text(size = 15),  axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
        facet_wrap(~sample_type, ncol = 5)

Subject as fixed effect vs strata term

Subject as fixed effect

adonis (dist ~ sample_type + log10(Final_reads) + treated + subject)

With strata

adonis (dist ~ sample_type + log10(Final_reads) + treated, strata = subject)

Strata term was employed rather than fixed effect

bray_perm_uni %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html", caption = "Subject id as fixed effect") %>%
        kable_styling(full_width = 0, html_font = "serif")

bray_perm_uni_strata %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html", caption = "Subject id as strata term") %>%
        kable_styling(full_width = 0, html_font = "serif")

Treatment significantly affected the beta-diversity.

Strata term is better representing our study aim.

What type of method affected the community at the most?

PERMANOVA on each treatment

dist ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp, strata = subject

QIAamp showed highest changes. But, it could be sample type specific.

bray_perm_strata %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

PERMANOVA with interaction term

dist ~ sample_type * treatment + log10(Final_reads), strata = subject

It was sample type specific. We need stratified analysis

bray_perm_inter %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "sample_type:treatment" ~ 'Sample type * treatment',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (Positive)

Untreateds chanced after treatment; by a lot. Differential abundance analysis should be conducted

bray_perm_pos %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (Neg)

Untreateds chanced after treatment

bray_perm_neg %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (Nasal)

lyPMA affected Nasal samples.

bray_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (BAL)

QIAamp affected BAL samples

bray_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (spt)

Sputum was affected by Molysis and QIAamp.

bray_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

MaAslin - volcano plot

Without interaction

feature ~ log10(final reads) + sample type + lyPMA + Benzonase + Host zero + Molysis + QIAamp + (1|subject)

Most samples are differentially abundant by sample type

#Making significance table for figure
        # Define a function to make species names italicized
species_italic <- function(data) {
  names <- gsub("_", " ", rownames(data))
  names <- gsub("[]]|[[]", "", names)
  names <- gsub(" sp", " sp.", names)
  names <- gsub(" sp.", "* sp.", names)
  names <- gsub(" group", "* group.", names)
  names <- ifelse(grepl("[*]", names), paste("*", names, sep = ""), paste("*", names, "*", sep = ""))
  rownames(data) <- names
  data
}

# Make a significance table for each figure (top 20 taxa)
make_sig_table <- function(data) {
  sig_data <- spread(data[order(data$qval), c("feature", "metadata", "qval")], metadata, qval)
  sig_data$min <- apply(sig_data, 1, FUN = min)
  sig_data <- sig_data[order(sig_data$min),] %>% select("feature", "lypma", "benzonase", "host_zero", "molysis", "qiaamp") %>% .[1:20,]
  sig_data[["feature"]] <- ifelse(sig_data[["feature"]] == "X.Collinsella._massiliensis", "[Collinsella]_massiliensis", sig_data[["feature"]])
  sig_data_italic <- sig_data %>% rownames_to_column(var = "-") %>%
          column_to_rownames(var = "feature") %>% species_italic %>% select(-c("-")) %>%
          rename(lyPMA = lypma,  Benzonase = benzonase, `Host zero` = host_zero, Molysis = molysis, QIAamp = qiaamp)
  sig_data_sig <- ifelse(sig_data_italic < 0.1, "*", NA) %>% data.frame(check.names = F)
  return(list(data = sig_data, data_italic = sig_data_italic, data_sig = sig_data_sig))
}

make_sig_table_neg <- function(data) {
  sig_data <- spread(fit_data_neg[order(fit_data_neg$qval), c("feature", "metadata", "qval")], metadata, qval)
  sig_data$min <- apply(sig_data, 1, FUN = min)
  sig_data <- sig_data[order(sig_data$min),] %>% select("feature", "lypma", "benzonase", "host_zero", "molysis", "qiaamp") %>% .[1:16,]
  sig_data[["feature"]] <- ifelse(sig_data[["feature"]] == "X.Collinsella._massiliensis", "[Collinsella]_massiliensis", sig_data[["feature"]])
  sig_data_italic <- sig_data %>% rownames_to_column(var = "-") %>%
          column_to_rownames(var = "feature") %>% species_italic %>% select(-c("-"))  %>% 
          rename(lyPMA = lypma,  Benzonase = benzonase, `Host zero` = host_zero, Molysis = molysis, QIAamp = qiaamp)
  sig_data_sig <- ifelse(sig_data_italic < 0.1, "*", NA) %>% data.frame(check.names = F)
  return(list(data = sig_data, data_italic = sig_data_italic, data_sig = sig_data_sig))
}

fit_data_neg <- make_sig_table_neg(fit_data_neg)
fit_data_pos <- make_sig_table(fit_data_pos)
fit_data_bal <- make_sig_table(fit_data_bal)
fit_data_ns <- make_sig_table(fit_data_ns)
fit_data_spt <- make_sig_table(fit_data_spt)
fit_data_bal_neg <- make_sig_table(fit_data_bal_neg)
fit_data_ns_neg <- make_sig_table(fit_data_ns_neg)
fit_data_spt_neg <- make_sig_table(fit_data_spt_neg)


neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Neg."),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Neg.")) %in% fit_data_neg$data$feature)
fit_data_neg$rel <- cbind(neg_sig %>% otu_table %>% t, neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_neg$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")


pos_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Mock"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Mock")) %in% fit_data_pos$data$feature)
fit_data_pos$rel <- cbind(pos_sig %>% otu_table %>% t, pos_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_pos$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")


spt_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Sputum")) %in% fit_data_spt$data$feature)
fit_data_spt$rel <- cbind(spt_sig %>% otu_table %>% t, spt_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_spt$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

spt_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "Sputum")) %in% fit_data_spt_neg$data$feature)

fit_data_spt_neg$rel <- cbind(spt_neg_sig %>% otu_table %>% t, spt_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_spt_neg$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")


ns_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Nasal"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Nasal")) %in% fit_data_ns$data$feature)
fit_data_ns$rel <- cbind(ns_sig %>% otu_table %>% t, ns_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>%
        .[row.names(fit_data_ns$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

ns_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Nasal"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "Nasal")
                                                  ) %in% fit_data_spt_neg$data$feature)
fit_data_ns_neg$rel <- cbind(ns_neg_sig %>% otu_table %>% t, ns_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_ns_neg$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")



bal_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "BAL"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "BAL")) %in% fit_data_bal$data$feature)
fit_data_bal$rel <- cbind(bal_sig %>% otu_table %>% t, bal_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_bal$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

bal_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "BAL"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "BAL")
                                                  ) %in% fit_data_bal_neg$data$feature)

fit_data_bal_neg$rel <- cbind(bal_neg_sig %>% otu_table %>% t, bal_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>% data.frame(check.names = F) %>% 
        .[row.names(fit_data_bal_neg$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

#Volcano plot

ggplot(maaslin_all, aes(y = -log10(qval), x = coef, col = metadata)) +
        theme_classic(base_family = "serif") +
        labs(tag = "A") +
        geom_point(size = 2) +
        xlab("MaAslin coefficient") +
        ylab("-log<sub>10</sub>(*q*-value)") +
        geom_hline(yintercept = 1, col = "gray") +
        geom_vline(xintercept = 0, col = "gray") +
        geom_richtext(aes( 4, 8, label = "*q*-value = 0.1, fold-change = 0", vjust = -1, fontface = 1), col = "grey", size = 3, family = "serif") +
        theme(legend.position = "top", axis.title.y = ggtext::element_markdown()) +
        scale_color_manual(values = c("#8c510a",  "#c51b7d", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"),
                           breaks = c("log10.Final_reads", "sample_type", "lypma", "benzonase", "host_zero",  "molysis", "qiaamp"), 
                           labels = c("log10 (Final reads)", "Sample type", "lyPMA", "Benzonase", "Host zero",  "Molysis", "QIAamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        guides(col = guide_legend(title = "Covariates", title.position = "top", nrow = 2))

MaAslin - table

feature ~ log10(final reads) + sample type + lyPMA + Benzonase + Host zero + Molysis + QIAamp + (1|subject)

Some taxa were changed due to treatment

Stratified analysis is required.

Table of associations had significant (q < 0.1) result

cat("all association")

## all association

maaslin_all %>% subset(., .$qval < 0.1) %>% arrange(., .$feature) %>% .$metadata %>% table

## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                27                15               100                25 
##           molysis            qiaamp       sample_type 
##                31                15               122

cat("Positive association (increased taxa)")

## Positive association (increased taxa)

maaslin_all %>% subset(., .$qval < 0.1) %>% subset(., .$coef > 0) %>% arrange(., .$feature) %>% .$metadata %>% table

## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                19                 7                99                24 
##           molysis            qiaamp       sample_type 
##                23                 7               111

MaAslin - can’t test global significance of a covariate with multi-level.

(https://forum.biobakery.org/t/global-significance-test-for-multilevel-factor/3061)

Baloon plot - Neg.

No taxa changed after treatment

#Mergeing data into one dataframe
merge(fit_data_neg$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_neg$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_neg$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Baloon plot - Pos.

Taxa decreased after treatment

#Mergeing data into one dataframe
merge(fit_data_pos$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_pos$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_pos$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Baloon plot - BAL

Mean relative abundances of top 20 taxa had low q-values.

No taxa changed after treatment

#Mergeing data into one dataframe
merge(fit_data_bal$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_bal$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_bal$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

#Mergeing data into one dataframe

merge(fit_data_bal_neg$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_bal_neg$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_bal_neg$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Baloon plot - Nasals

Mean relative abundances of top 20 taxa had low q-values.

Some taxa changed after treatment, but nothing was unique

#Mergeing data into one dataframe
merge(fit_data_ns$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_ns$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_ns$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

#Data with negative
merge(fit_data_ns_neg$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_ns_neg$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_ns_neg$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Baloon plot - Sputum

Mean relative abundances of top 20 taxa had low q-values.

Some taxa changed after treatment, but nothing was unique

#Mergeing data into one dataframe
merge(fit_data_spt$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_spt$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_spt$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

#Mergeing data into one dataframe

merge(fit_data_spt_neg$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      fit_data_spt_neg$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(fit_data_spt_neg$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Results

Some taxa were significantly changed after treatment. Among top 20, no taxa observed in only one treatment group. As their emergence were consistent across all treatment groups, they were considered as endogenus.

MaAslin with interaction

feature ~ log10(Final reads) + treatment + sample type + treatment * sample type (1|subject)

Some taxa were treaetment specific, after adjusting interaction of sample type * treatment

#Generating interaction term
#sample_data(phyloseq_rel_nz)$sampletype_treatment <- paste(sample_data(phyloseq_rel_nz)$sample_type, #sample_data(phyloseq_rel_nz)$treatment, sep = "*")

#capture.output(maaslin_interaction = Maaslin2(input_data = otu_table(phyloseq_rel_nz) %>% t %>% data.frame(), 
#                 input_metadata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), 
#                 output = "data", 
#                 fixed_effects = c("sample_type", "log10.Final_reads", "treatment", "sampletype_treatment"), 
#                 transform = "LOG", #default
#                 normalization = "TSS", 
#                 random_effects = c("subject_id"), 
#                 reference = c("sample_type,BAL", "treatment,Untreated", "sampletype_treatment,BAL*Untreated"), 
#                 plot_heatmap = F,
#                 plot_scatter = F))

maaslin_interaction <- read.csv("data/maaslin_interaction.csv")

#interaction term - ggplot
ggplot(maaslin_interaction, aes(y = -log10(qval), x = coef, col = metadata)) +
         theme_classic(base_family = "serif") +
         #labs(tag = "A") +
         ggtitle("MaAslin with interaction term")+
         geom_point(size = 2) +
         xlab("MaAslin coefficient") +
         ylab("-log<sub>10</sub>(*q*-value)") +
         geom_hline(yintercept = 1, col = "gray") +
         geom_vline(xintercept = 0, col = "gray") +
         geom_richtext(aes( 4, 8, label = "*q*-value = 0.1, fold-change = 0", vjust = -1, fontface = 1), col = "grey", size = 3, family = "serif") +
         theme(legend.position = "top", axis.title.y = ggtext::element_markdown()) +
         scale_color_manual(values = c("#e41a1c",  "#377eb8", "#4daf4a", "#984ea3")) +
         guides(col = guide_legend(title = "Fixed effects", title.position = "top", nrow = 1))

#Checking number of bugs differentially abundance with interaction term 
cat("Number of differentially abundant bugs by each metadata")

## Number of differentially abundant bugs by each metadata

maaslin_interaction %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()

## .
##    log10.Final_reads          sample_type sampletype_treatment 
##                   38                  146                  560 
##            treatment 
##                  135

MaAsLin interaction analysis

Hypothesis: if a sample is contaminated by some treatment, a change of taxon is likely to be associated with one treatment method

No taxa increased only because of one treatment method

 cat("Some taxa were increased by each treatmment.\n But they are not contaminants, \nif they are present in most of the treatments")

## Some taxa were increased by each treatmment.
##  But they are not contaminants, 
## if they are present in most of the treatments

 maaslin_interaction %>% subset(., .$qval < 0.1 & .$metadata == "treatment") %>% .$feature %>% table %>% data.frame %>% arrange(-Freq) %>% rename(Feature = ".") %>% kbl(format = "html", caption = "Table of taxa differentially abundant by treatment") %>%
        kable_styling(full_width = 0, html_font = "serif")

 cat("Most of taxa were found on most of treatments.")

## Most of taxa were found on most of treatments.

 cat("Some taxa were treatment specific, only to one group")

## Some taxa were treatment specific, only to one group

subset(maaslin_interaction, maaslin_interaction$feature %in%  (maaslin_interaction %>% subset(., .$qval < 0.1 & .$metadata == "treatment") %>%
         .$feature %>% table %>% data.frame %>% subset(., Freq == 1) %>% .$. %>% as.character())) %>% subset(., .$qval < 0.1) %>% select(c("feature", "metadata", "value", "coef", "qval")) %>% remove_rownames() %>% kbl(format = "html", caption = "Table of taxa specific to one treatment group") %>%
        kable_styling(full_width = 0, html_font = "serif")

No taxa was increased due to one treatmemnt.

decontam - all sample

Listeria floridensis could be a potential contaminant

# Decontam package --------------------------------------------------------

# common contaminants across all the treatment methods
#Decontam - were there any contaminants?#
sample_data(phyloseq$phyloseq_rel)$is.neg <- grepl("Neg", sample_data(phyloseq$phyloseq_rel)$sample_type)
phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_rel, S.obs != 0)
#With all sampels
dec_f_all <- isContaminant(phyloseq_rel_nz, method="frequency", conc="DNA_bac_well")
dec_p_all <- isContaminant(phyloseq_rel_nz, method="prevalence", neg="is.neg", threshold=0.5)
dec_c_all <- isContaminant(phyloseq_rel_nz, method="combined", neg="is.neg", conc = "DNA_bac_well")

cat("decontam frequency - all sample")

## decontam frequency - all sample

dec_f_all %>% subset(.,.$contaminant)

cat("decontam prevalence - all sample")

## decontam prevalence - all sample

dec_p_all %>% subset(.,.$contaminant)

cat("decontam combined - all sample")

## decontam combined - all sample

dec_c_all %>% subset(.,.$contaminant)

decontam - stratified by sample_type

Stratified analysis showed no contaminants in NS and BAL

Sputum may have Corynebacterium pseudodiphtheriticum and Candida albicans as contaminants.

#Stratified by sample type

cat("decontam prevalence - BAL")

## decontam prevalence - BAL

subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)

cat("decontam prevalence - Nasal")

## decontam prevalence - Nasal

subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal", "Neg")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)

cat("decontam prevalence - Sputum")

## decontam prevalence - Sputum

subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)

cat("decontam frequency - BAL")

## decontam frequency - BAL

subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam frequency - Nasal")

## decontam frequency - Nasal

subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal", "Neg")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam frequency - Sputum")

## decontam frequency - Sputum

subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - BAL")

## decontam combined - BAL

subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Nasal")

## decontam combined - Nasal

subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal", "Neg")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Sputum")

## decontam combined - Sputum

subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

Stratified analysis showed no contaminants in NS and BAL

Sputum may have Corynebacterium pseudodiphtheriticum and Candida albicans as contaminants.

#Stratified by treatment

cat("decontam prevalence - lyPMA")

## decontam prevalence - lyPMA

subset_samples(phyloseq_rel_nz, treatment = "lypma") %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)

cat("decontam prevalence - lyPMA")

## decontam prevalence - lyPMA

subset_samples(phyloseq_rel_nz, treatment = "lypma") %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)

cat("decontam prevalence - lyPMA")

## decontam prevalence - lyPMA

subset_samples(phyloseq_rel_nz, treatment = "lypma") %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)

cat("decontam frequency - Benzonase")

## decontam frequency - Benzonase

subset_samples(phyloseq_rel_nz, treatment = "benzonase") %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam frequency - Benzonase")

## decontam frequency - Benzonase

subset_samples(phyloseq_rel_nz, treatment = "benzonase") %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam frequency - Benzonase")

## decontam frequency - Benzonase

subset_samples(phyloseq_rel_nz, treatment = "benzonase") %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Molysis")

## decontam combined - Molysis

subset_samples(phyloseq_rel_nz, treatment = "molysis") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Molysis")

## decontam combined - Molysis

subset_samples(phyloseq_rel_nz, treatment = "molysis") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Molysis")

## decontam combined - Molysis

subset_samples(phyloseq_rel_nz, treatment = "molysis") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Host zero")

## decontam combined - Host zero

subset_samples(phyloseq_rel_nz, treatment = "host_zero") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Host zero")

## decontam combined - Host zero

subset_samples(phyloseq_rel_nz, treatment = "host_zero") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - Host zero")

## decontam combined - Host zero

subset_samples(phyloseq_rel_nz, treatment = "host_zero") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - QIAamp")

## decontam combined - QIAamp

subset_samples(phyloseq_rel_nz, treatment = "qiaamp") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - QIAamp")

## decontam combined - QIAamp

subset_samples(phyloseq_rel_nz, treatment = "qiaamp") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

cat("decontam combined - QIAamp")

## decontam combined - QIAamp

subset_samples(phyloseq_rel_nz, treatment = "qiaamp") %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

Function richness (all samples & interaction term) - ANOVA

Interaction term showed high p values. However, it could be due to even effect sample type * treatment. Interaction term will be tested.

sample_data <- sample_data(phyloseq$phyloseq_path_rpkm) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))
lmer_sob <- lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data)
lmer_sob %>% 
        anova() %>% 
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Function richness (all samples & interaction term)

Effect of some treatment was neutralized by interactin term. Therefore, the association was sample_type specific.

Stratified analysis will be conducted.

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

All terms were significant.

Function richness - stratified (NS)

Some treatment enabled discovering more functions in Nasals

lmer(S.obs ~ treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Function richness (BAL)

Higher Final reads enables more discovery of functions.

lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Function richness (sputum)

Sputum showed no changes. This may due to an enrichment of richness in control groups.

lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inv Simp - ANOVA

p - value = 0.096 for the interaction term. Interaction term will be tested.

lmer_invsimpson <- lmer(data_invsimpson ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data)
lmer_invsimpson %>% 
        anova() %>% 
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inv. Simpson (all samples & interaction term)

Sample type specific effect was observed. Stratified anlysis required.

#Simpson

lmer_invsimpson %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (NS)

Inverse Simpson ~ sample_type + log10 (Final_reads) + (1|original_sample)

Nasal showed changes after Molysis treatment

lmer(data_invsimpson ~ treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (BAL)

No changes found at BAL

lmer(data_invsimpson ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (spt)

Changes associated with deeper sequencing with sputum

lmer(data_invsimpson ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Treatment as categorized group

dist ~ sample_type + log10(Final_reads) + treatment, strata = subject

No significant changes were observed.

bray_perm_uni_strata %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html", caption = "Subject id as strata term") %>%
        kable_styling(full_width = 0, html_font = "serif")

Function PERMAONVA (detailed treatment)

dist ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp, strata = subject

No significant changes were observed.

bray_perm_strata %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

QIAamp showed highest changes. But, it could be sample type specific.

Function interaction term

dist ~ sample_type * treatment + log10(Final_reads), strata = subject

Some changes were treatment induced?

We don’t have to run stratified anlysis

bray_perm_inter %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "sample_type:treatment" ~ 'Sample type * treatment',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (NS)

Stratified analysis not matching with boxplot results

bray_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (BAL)

Stratified analysis not matching with boxplot results

bray_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Stratified (spt)

Stratified analysis not matching with boxplot results

bray_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

MaAslin without interaction - volcano plot

Again, most of DA functions were sample type specific

#Making significance table for figure
        # Define a function to make species names italicized
# Make a significance table for each figure (top 20 taxa)
make_sig_table <- function(data) {
  sig_data <- spread(data[order(data$qval), c("feature", "metadata", "qval")], metadata, qval)
  sig_data$feature <- gsub("[.]", "-", sig_data$feature)
  sig_data$min <- apply(sig_data, 1, FUN = min)
  sig_data <- sig_data[order(sig_data$min),] %>% select("feature", "lypma", "benzonase", "host_zero", "molysis", "qiaamp") %>% .[1:20,]
  sig_data[["feature"]] <- ifelse(sig_data[["feature"]] == "X.Collinsella._massiliensis", "[Collinsella]_massiliensis", sig_data[["feature"]])
  sig_data_italic <- sig_data %>% rownames_to_column(var = "-") %>%
          column_to_rownames(var = "feature") %>% select(-c("-")) %>%
          rename(lyPMA = lypma,  Benzonase = benzonase, `Host zero` = host_zero, Molysis = molysis, QIAamp = qiaamp)
  sig_data_sig <- ifelse(sig_data_italic < 0.1, "*", NA) %>% data.frame(check.names = F)
  return(list(data = sig_data, data_italic = sig_data_italic, data_sig = sig_data_sig))
}

f_fit_data_neg <- make_sig_table(f_fit_data_neg)
f_fit_data_pos <- make_sig_table(f_fit_data_pos)
f_fit_data_bal <- make_sig_table(f_fit_data_bal)
f_fit_data_ns <- make_sig_table(f_fit_data_ns)
f_fit_data_spt <- make_sig_table(f_fit_data_spt)
f_fit_data_bal_neg <- make_sig_table(f_fit_data_bal_neg)
f_fit_data_ns_neg <- make_sig_table(f_fit_data_ns_neg)
f_fit_data_spt_neg <- make_sig_table(f_fit_data_spt_neg)



f_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Neg."),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Neg.")) %in% f_fit_data_neg$data$feature)
f_fit_data_neg$rel <- cbind(f_neg_sig %>% otu_table %>% t, f_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>% 
        .[row.names(f_fit_data_neg$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")


f_pos_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Mock"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Mock")) %in% f_fit_data_pos$data$feature)
f_fit_data_pos$rel <- cbind(f_pos_sig %>% otu_table %>% t, f_pos_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>% 
        .[row.names(f_fit_data_pos$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")



f_spt_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Sputum")) %in% f_fit_data_spt$data$feature)
fit_data_spt$rel <- cbind(f_spt_sig %>% otu_table %>% t, f_spt_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>% 
        .[row.names(f_fit_data_spt$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

f_spt_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "Sputum")) %in% f_fit_data_spt$data$feature)

f_fit_data_spt$rel <- cbind(f_spt_sig %>% otu_table %>% t, f_spt_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>% 
        .[row.names(f_fit_data_spt$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

f_spt_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "Sputum")) %in% f_fit_data_spt_neg$data$feature)

f_fit_data_spt_neg$rel <- cbind(f_spt_neg_sig %>% otu_table %>% t, f_spt_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>% 
        .[row.names(f_fit_data_spt_neg$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")


f_ns_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Nasal"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "Nasal")) %in% f_fit_data_ns$data$feature)

f_fit_data_ns$rel <- cbind(f_ns_sig %>% otu_table %>% t, f_ns_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>% 
        .[row.names(f_fit_data_ns$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")
f_fit_data_ns$rel$feature <- row.names(f_fit_data_ns$data_sig)


f_ns_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "Nasal"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "Nasal")) %in% f_fit_data_spt_neg$data$feature)

f_fit_data_ns_neg$rel <- cbind(f_ns_neg_sig %>% otu_table %>% t, f_ns_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>%
        .[row.names(f_fit_data_spt_neg$data_italic),] %>%
        mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")



f_bal_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "BAL"),
                                       taxa_names(subset_samples(phyloseq_rel_nz, sample_type == "BAL")) %in% f_fit_data_bal$data$feature)

f_fit_data_bal$rel <- cbind(f_bal_sig %>% otu_table %>% t, f_bal_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>%
        .[row.names(f_fit_data_bal$data_italic),] %>%
        mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")




f_bal_neg_sig <- subset_taxa(subset_samples(phyloseq_rel_nz, sample_type == "BAL"),
                                       taxa_names(subset_samples(phyloseq_rel_nz,
                                                                 sample_type == "BAL")
                                                  ) %in% f_fit_data_bal_neg$data$feature)

f_fit_data_bal_neg$rel <- cbind(f_bal_neg_sig %>% otu_table %>% t, f_bal_neg_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% data.frame(check.names = F) %>%
        .[row.names(f_fit_data_bal_neg$data_italic),] %>%
        mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")




#Volcano plot

ggplot(f_maaslin_all, aes(y = -log10(qval), x = coef, col = metadata)) +
        theme_classic(base_family = "serif") +
        labs(tag = "A") +
        geom_point(size = 2) +
        xlab("MaAslin coefficient") +
        ylab("-log<sub>10</sub>(*q*-value)") +
        geom_hline(yintercept = 1, col = "gray") +
        geom_vline(xintercept = 0, col = "gray") +
        #geom_richtext(aes( 4, 8, label = "*q*-value = 0.1, fold-change = 0", vjust = -1, fontface = 1), col = "grey", size = 3, family = "serif") +
        theme(legend.position = "top", axis.title.y = ggtext::element_markdown()) +
        scale_color_manual(values = c("#4daf4a",  "#984ea3", "#f781bf", "#377eb8", "#ff7f00", "#ffff33", "#a65628"),
                           breaks = c("log10.Final_reads", "sample_type", "lypma", "benzonase", "host_zero",  "molysis", "qiaamp"), 
                           labels = c("log10 (Final reads)", "Sample type", "lyPMA", "Benzonase", "Host zero",  "Molysis", "QIAamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        guides(col = guide_legend(title = "Covariates", title.position = "top", nrow = 2))

Most of the DA function were sample type dependent.

MaAsLin table (function)

Large number of functions were differentially aubundant.

f_maaslin_all %>% subset(., .$qval < 0.1) %>% .$metadata %>% table

## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                90               100               368                83 
##           molysis            qiaamp       sample_type 
##                16                36               591

Stratified analysis is required.

Baloon plot - Nasals

Similarly, few functions were newly discovered

merge(f_fit_data_spt$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_spt$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(f_fit_data_spt$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Baloon plot - BAL

Similarly, few functions were newly discovered

merge(f_fit_data_bal$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_bal$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(f_fit_data_bal$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Baloon plot - Nasal

Some functions were newly discovered.

merge(f_fit_data_ns$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_ns$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(f_fit_data_ns$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
#Baloon plot
        ggballoonplot(size = "value", fill = "qval", y = "feature", x= "treatment") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(vjust = 0.5, hjust=0.5),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown())  +
        #Adding significance asterisks
        geom_text(aes(y = feature,
                      x = treatment,
                      label = sig,
                      col = "red"),
                  hjust = -2,
                  vjust = 0.8,
                  size = 5) +
        guides(col = guide_legend(nrow = 1,
                                  override.aes = aes(label = "*", size = 10, color = "red"),
                                  title="Significance",
                                  title.position = "top", order = 3, ),
               fill = guide_colorbar(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top",
                                   order = 2),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 2),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         ) +
        scale_color_manual(values = c("red"),
                           labels = c(expression(paste(italic("q"),
                                                       "-value < 0.1",
                                                       sep = "")
                                                 )
                                      )
                           )

Results After adding control data, MaAslin needs to be reanalyzed. Adding controls (mock communities) for each treatment group will show more statistically valid results in y ~ log(final reads) + sample_type + treatment, (re = subject_id))